function testoutputdir_full=osl_test_script_group_full(testdatadir,testoutputdir,S)

global OSLDIR;

hmm_av_class_occupancy=20;
hmm_do_glm_statewise=1;
forward_meg=S.forward_meg;
do_hmm=S.do_hmm;
do_tf=S.do_tf;

testoutputdir_full=[testoutputdir '/script_group_full_hmmav' num2str(hmm_av_class_occupancy) '_' regexprep(S.forward_meg, ' ', '') '_hmm' num2str(S.do_hmm) '_' S.recon_method '_sss' num2str(S.do_sss) '_tf' num2str(S.do_tf)];
testoutputdir_full

runcmd(['rm -rf ' testoutputdir_full]);
mkdir(testoutputdir_full);

testplotsdir=[testoutputdir_full '/plots'];
mkdir(testplotsdir);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% osl_example_group_oat

printprefix='group_oat';
printindex=1;

datadir=[testdatadir '/faces_group_data_new_full']; % this is the directory the oat files will be stored in
%if(~S.do_sss)
%    datadir=['/home/disk3/mwoolrich/vols_data/murphy_project/data']; % this is the directory the oat files will be stored in
%else
%    datadir=['/home/disk3/mwoolrich/vols_data/osl_testdata_dir/faces_group_data_new_full'];
%end;

spmfilesdir=[datadir '/spm_files']; % this is the directory the SPM files will be stored in
structuralsdir=[datadir '/structurals']; % this is the directory the structural files will be found in

clear fif_files spm_files structural_files

% Alternatively you can leave fif_files empty and set up a list of SPM MEEG
% object files:
is=[15:28];
spm_files=[];
for i=1:length(is),  
    if(S.do_sss)
        spm_files{i}=[spmfilesdir '/espm8_meg' num2str(is(i)) '.mat'];
    else
        spm_files{i}=[spmfilesdir '/efspm8_nosss_meg' num2str(is(i)) '_1.mat'];
    end;
end;

% structural files:
sfs={'F9'    'F11'    'F15'    'F16'    'F20'    'F22'    'F25'    'M3'    'M6'    'M7'    'M12'    'M14'    'M18'    'M19'    'F4'    'F7'    'F10'    'F13'    'F18'    'F24'    'F26'    'M2'    'M4'    'M8'    'M10'    'M13'   'M17'    'M21'    'F3'    'F6'    'F12'    'F14'    'F17'    'F19'    'F23'    'M5'    'M9'    'M11'    'M15'    'M20'    'M22'    'M23'};
sfs{21}=''; % structural is missing
sfs{14}=''; % structural is missing
sfs{6}=''; % structural is missing
for i=1:length(is),
    if ~strcmp(sfs{is(i)},'')
        structural_files{i}=[structuralsdir '/' sfs{is(i)} '/struct.nii'];
    else
        structural_files{i}='';
    end;
end;

oat=[];
oat.source_recon.D_continuous=[]; % do not have continuous files
oat.source_recon.D_epoched=spm_files; % only epoched files
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
if do_tf,
    oat.source_recon.freq_range=[35 80]; % frequency range in Hz
    oat.source_recon.time_range=[-0.2 0.6];
    oat.first_level.time_range=[-0.1 0.5];
    oat.first_level.tf_num_freqs=1;
    oat.first_level.tf_method='hilbert';
else
    oat.source_recon.freq_range=[1 48]; % frequency range in Hz
    oat.source_recon.time_range=[-0.2 0.35];
    oat.first_level.num_freqs=1;
    oat.first_level.tf_method='none';
end;

oat.source_recon.method=S.recon_method;
oat.source_recon.gridstep=7; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.mri=structural_files;
oat.source_recon.sessions_to_do=1:14;
oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat.source_recon.forward_meg=forward_meg;
oat.source_recon.pca_dim=58;
oat.source_recon.force_pca_dim=1;

if(do_hmm)
    oat.source_recon.hmm_num_states=-1;
    oat.source_recon.hmm_num_starts=4;
    oat.source_recon.hmm_pca_dim=40;
    oat.source_recon.hmm_av_class_occupancy=hmm_av_class_occupancy;
end;

Xsummary={};Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast{4}=[0 0 -1 1]'; 
oat.first_level.contrast{5}=[0 -1 0 1]';
oat.first_level.contrast{6}=[3 1 1 1]';
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
oat.first_level.contrast_name{4}='fear-happy';
oat.first_level.contrast_name{5}='fear-neutral';
oat.first_level.contrast_name{6}='faces+motorbikes';
oat.first_level.cope_type='acope';
oat.first_level.hmm_do_glm_statewise=hmm_do_glm_statewise;

oat = osl_check_oat(oat);

% group-level stage

%oat = osl_load_oat(oat.source_recon.dirname,'first_level','sub_level','group_level'); 
oat.source_recon.sessions_to_do=1:length(spm_files);
oat.first_level.sessions_to_do=oat.source_recon.sessions_to_do;
oat.subject_level.subjects_to_do=oat.first_level.sessions_to_do;

%oat.group_level.time_range=[0 0.3];
oat.group_level.space_average=0;
oat.group_level.time_average=0;
oat.group_level.time_smooth_std=0; % secs
oat.group_level.spatial_smooth_fwhm=0; % mm
oat.group_level.mask_fname='';oat.group_level=rmfield(oat.group_level,'mask_fname');
oat.group_level.use_tstat=0;
oat.group_level.group_varcope_time_smooth_std=2;
oat.group_level.group_varcope_spatial_smooth_fwhm=10;
oat.group_level.name='group_level';

% run OAT
oat.to_do=[1 1 1 1]; % run group-level stage only

oat = osl_check_oat(oat);

oat = osl_run_oat(oat);

if(0)
    % OUTPUT SUBJECT 1'S NIFTII FILES

    str=['fslview ' OSLDIR '/std_masks/MNI152_T1_' num2str(oat.source_recon.gridstep) 'mm_brain '];
    for ii=1:1,%14,
        S2=[];
        S2.oat=oat;
        S2.stats_fname=oat.first_level.results_fnames{ii};
        S2.first_level_contrasts=[3]; % list of first level contrasts to output
        S2.resamp_gridstep=oat.source_recon.gridstep;

        [statsdir,times]=osl_save_nii_stats(S2);
        str=[str statsdir '/tstat3_10mm '];
    end;

    str2=[str ' &'];
    runcmd(str2);
   
end;


%% OUTPUT GROUP'S NIFTII FILES

S2=[];
S2.oat=oat;
S2.stats_fname=oat.group_level.results_fnames;
S2.first_level_contrasts=[1:6]; % list of first level contrasts to output

[statsdir,times]=osl_save_nii_stats(S2);

figure;
con=3;
map=ra([statsdir '/tstat' num2str(con) '_gc1_2mm']);
resamp_gridstep=2;bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=99.5;percto=99.9;
low=percentile((x1),percfrom);high=percentile((x1),percto);
zslice=30;vol=44;
overlay_act(flipud(squeeze(map(:,:,zslice,vol))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
text(2,2,[num2str(low) ',' num2str(high)],'Color','w');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

% !cp /Users/woolrich/homedir/matlab/osl_testdata_dir/test_output_script_group_full_hmmav20_MEGLocalSpheres_hmm0/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /Users/woolrich/homedir/matlab/osl_testdata_dir/tstat_full_localsphere_hmm0.nii.gz
% !cp /Users/woolrich/homedir/matlab/osl_testdata_dir/test_output_script_group_full_hmmav20_MEGLocalSpheres_hmm1/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /Users/woolrich/homedir/matlab/osl_testdata_dir/tstat_full_localsphere_hmm1.nii.gz
% !cp /Users/woolrich/homedir/matlab/osl_testdata_dir/test_output_script_group_full_hmmav20_SingleShell_hmm0/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /Users/woolrich/homedir/matlab/osl_testdata_dir/tstat_full_shell_hmm0.nii.gz
% !cp /Users/woolrich/homedir/matlab/osl_testdata_dir/test_output_script_group_full_hmmav20_SingleShell_hmm1/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /Users/woolrich/homedir/matlab/osl_testdata_dir/tstat_full_shell_hmm1.nii.gz
% !fslview /Users/woolrich/homedir/matlab/osl_testdata_dir/tstat_full_shell_hmm1.nii.gz /Users/woolrich/homedir/matlab/osl_testdata_dir/tstat_full_shell_hmm0.nii.gz /Users/woolrich/homedir/matlab/osl_testdata_dir/tstat_full_localsphere_hmm1.nii.gz /Users/woolrich/homedir/matlab/osl_testdata_dir/tstat_full_localsphere_hmm0.nii.gz


% !cp /home/disk3/mwoolrich/vols_data/osl_testdata_dir/test_output_script_group_full_hmmav20_MEGLocalSpheres_hmm0/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_localsphere_hmm0.nii.gz
% !cp /home/disk3/mwoolrich/vols_data/osl_testdata_dir/test_output_script_group_full_hmmav20_MEGLocalSpheres_hmm1/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_localsphere_hmm1.nii.gz
% !cp /home/disk3/mwoolrich/vols_data/osl_testdata_dir/test_output_script_group_full_hmmav20_SingleShell_hmm0/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_shell_hmm0.nii.gz
% !cp /home/disk3/mwoolrich/vols_data/osl_testdata_dir/test_output_script_group_full_hmmav20_SingleShell_hmm1/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_shell_hmm1.nii.gz
% !fslview /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_shell_hmm1.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_shell_hmm0.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_localsphere_hmm1.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_localsphere_hmm0.nii.gz

% !cp /home/disk3/mwoolrich/vols_data/osl_testdata_dir/test_output_script_group_full_MEGLocalSpheres_hmm0/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_localsphere_hmm0.nii.gz
% !cp /home/disk3/mwoolrich/vols_data/osl_testdata_dir/test_output_script_group_full_MEGLocalSpheres_hmm1/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_localsphere_hmm1.nii.gz
% !cp /home/disk3/mwoolrich/vols_data/osl_testdata_dir/test_output_script_group_full_SingleShell_hmm0/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_shell_hmm0.nii.gz
% !cp /home/disk3/mwoolrich/vols_data/osl_testdata_dir/test_output_script_group_full_SingleShell_hmm1/group_oat.oat/first_level_sub_level_group_level_dir/tstat3_gc1_2mm.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_shell_hmm1.nii.gz
% !fslview /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_shell_hmm1.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_shell_hmm0.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_localsphere_hmm1.nii.gz /home/disk3/mwoolrich/vols_data/osl_testdata_dir/tstat_full_localsphere_hmm0.nii.gz


for cons=1:length(S2.first_level_contrasts),
    con=S2.first_level_contrasts(cons);
    runcmd(['cp ' testoutputdir_full '/group_oat.oat/first_level_sub_level_group_level_dir/tstat' num2str(con) '_gc1_2mm.nii.gz ' testoutputdir '/tstat' num2str(con) '_full_' regexprep(S.forward_meg, ' ', '') '_hmm' num2str(S.do_hmm) '_' S.recon_method '_sss' num2str(S.do_sss) '.nii.gz']);
end;
    
%%%%%%%%%%%%%%%%%%
%% Using ROI at the input to the group level

% load OAT analysis for which the first 3 stages have already been run
oatdir=oat.source_recon.dirname;

oat = osl_load_oat(oatdir,'first_level','sub_level','group_level'); 

% mask
oat.group_level.time_range=[-0.2 0.3];
oat.group_level.space_average=1;
oat.group_level.time_average=0;
oat.group_level.time_smooth_std=0; % secs
%oat.group_level.group_varcope_time_smooth_std=0.2;
%oat.group_level.group_varcope_spatial_smooth_fwhm=100;
oat.group_level.spatial_smooth_fwhm=0; % mm
oat.group_level.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
oat.group_level.mask_fname=[OSLDIR '/std_masks/Right_Amygdala'];

oat.group_level.store_lower_level_copes=1;

oat.group_level.name='roi_average_group_level';
oat.group_level.use_tstat=0;

oat = osl_check_oat(oat);

oat.to_do=[0 0 0 1];

% run OAT
oat = osl_run_oat(oat);

%% do plots

gstats=osl_load_oat_results(oat,oat.group_level.results_fnames);

tim=gstats.times
cons=[1:5];

figure;
baseline_correct=0;
gcon=1;
subplot(1,2,1);ho;
cols={'r','g','b','k','m'};
for coni=1:length(cons),
    con=cons(coni);
    cope=(permute(mean(gstats.cope(:,:,con,:,gcon),1),[2 4 1 3 5]));
    stdcope=(permute(mean(gstats.stdcope(:,:,con,:,gcon),1),[2 4 1 3 5]));
    mn=0;if(baseline_correct),mn=mean(cope(tim>-0.1 & tim<0));end;
    errorbar(tim,cope-mn,stdcope,cols{coni});
end;
plot4paper('time (secs)','COPE'); 
title(['Group average copes (ERFs), use_tstat=' num2str(oat.group_level.use_tstat)]);
legend(oat.first_level.contrast_name,'NorthWest'); 
%exportfig(gcf,[murphy_plots_dir '/' Sin.resdir_postfix '_con' num2str(Sin.contrast4localisation) '_' Sin.mask_fname '/p_fan_cope'],'Color','rgb','Format','eps');
 
% tstat
baseline_correct=0;
gcon=1;
subplot(1,2,2);ho;
for coni=1:length(cons),
    con=cons(coni);
    cope=permute(mean(gstats.cope(1,:,con,:,gcon),1),[2 4 1 3 5]);
    stdcope=permute(mean(gstats.stdcope(1,:,con,:,gcon),1),[2 4 1 3 5]);
    mn=0;if(baseline_correct),mn=mean(cope(tim>-0.1 & tim<0));end;
    plot(tim,(cope-mn)./stdcope,cols{coni});
end;
plot4paper('time (secs)','t-stat'); 
title(['Group average t-stats (ERFs), use_tstat=' num2str(oat.group_level.use_tstat)]);
legend(oat.first_level.contrast_name,'NorthWest'); 

%% individual subject plots:

baseline_correct=1;
figure;
cs=[1:5];
for c=1:length(cs),
    subplot(2,ceil(length(cs)/2),c)
    con=cs(c);
    tmp=((gstats.lower_level_copes{con}));
    mn=0;if(baseline_correct),mn=mean(squash(tmp(1,:,tim<0)));end;
    plot(tim,permute(tmp(1,:,:),[2,3,1])-mn);ho;
    plot(tim,squeeze(median(tmp(1,:,:)-mn,2)),'LineWidth',3);
    plot4paper('time (secs)','1st level cope'); 
    title(['1st lev cope ' num2str(con)]);
end;

% figure;plot(tim,squeeze(mean(tmp(1,:,:),2))./sqrt(squeeze(var(tmp(1,:,:),[],2)./(size(tmp,2)-1)))); title('tstat');

%%%%%%%%%%%%%%%%%%%%%%%%
%% test session averaging

if(0),
    oat = osl_load_oat(oat.source_recon.dirname,'first_level','sub_level','group_level'); 
    oat.subject_level.session_index_list=[];
    oat.subject_level.session_index_list{1}=1:7;
    oat.subject_level.session_index_list{2}=8:14;
    oat.subject_level.subjects_to_do=[1 2];
    oat.subject_level.name='sub_name_session_av';

    oat.to_do=[0 0 1 0]; % run group-level stage only

    oat = osl_run_oat(oat);

    S2=[];
    S2.oat=oat;
    S2.stats_fname=oat.subject_level.results_fnames{1};
    S2.first_level_contrasts=[3]; % list of first level contrasts to output
    [statsdir,times]=osl_save_nii_stats(S2);

    figure;
    con=3;
    map=ra([statsdir '/tstat' num2str(con) '_2mm']);
    resamp_gridstep=2;bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
    x1=squash(map,abs(map));
    percfrom=99;percto=99.9;
    low=percentile((x1),percfrom);high=percentile((x1),percto);
    zslice=30;vol=44;
    overlay_act(flipud(squeeze(map(:,:,zslice,vol))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
    text(2,2,[num2str(low) ',' num2str(high)],'Color','w');
    print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;
end;

disp('*************************************************');
disp('Finished osl_test_scrip_group_full test');
disp('*************************************************');